Forcing reversibility in the no strand-bias substitution model allows for the theoretical 
and practical identifiability of its 5 parameters from pairwise DNA sequence 

comparisons. 

Osvaldo ZagordQ 
International School of Advanced Studies SISSA-ISAS 
via Beirut 2-4, 34013 Trieste, Italy 

Jean R. Lobry 
Laboratoire BBE-CNRS-UMR-5558, Univ. C. Bernard - Lyon I 
43 Bd 11/11/1918, F-69622 Villeurbanne CEDEX, France 

Because of the base pairing rules in DNA, some mutations experienced by a portion of DNA 
during its evolution result in the same substitution, as we can only observe differences in coupled 
nucleotides. Then, in the absence of a bias between the two DNA strands, a model with at most 
6 different parameters instead of 12 is sufficient to study the evolutionary relationship between 
homologous sequences derived from a common ancestor. On the other hand the same symmetry 
reduces the number of independent observations which can be made. Such a reduction can in some 
cases invalidate the calculation of the parameters. A compromise between biologically acceptable 
hypotheses and tractability is introduced and a five parameter reversible no- strand-bias condition 
(RNSB) is presented. The identifiability of the parameters under this model is shown by examples. 



I. INTRODUCTION 



Darwinian Evolution is based upon the interplay of two driving forces: mutation of an organism features, and 
natural selection acting on the living organisms. Nowadays the role of the DNA in the evolutive processes has been 
recognised, and the physical basis of the mutation process has been identified. Mutation acts on the DNA and we 
call mutation rate the probability that a descendant has a difference in the genome if this is compared to that of its 
parents. The substitution rate is the probability of finding a difference when comparing the genomes of species to one 
of its ancestors. 

We see that while the mutation is closely related to the biophysical process of DNA damage, or replication error 
etc., the substitution is the result of a mutation and of a population-dynamics process, which has spread the former 
to the whole population. A fundamental observation by M. Kimura in 1968 [y| argued that, in the case of neutral 
mutations (i.e. those mutations which have no apparent effect on the adaptation of an organism to the environment), 
we can deduce the mutation rate from the substitutions, as they are actually the same. 

Let's consider an ancestor O at time i = which separates into two different evolutive lineages, resulting in two 
different species, A and B at time t. It would be useful to define a distance between A and B and to have a tool to 
calculate it by just comparing the genomes of A and B. 

In order to study evolutionary distance between homologous DNA sequences (descending from a common ancestor) 
and their consequent relationship, a model for nucleotide substitution can be introduced. Generally, the process is 
assumed to be a Markov chain, if some assumptions are made about the underlying process. The general hypotheses 
are: 

• substitution rates do not depend on the position along the DNA sequence; 

• they are constant during evolutionary time; 

• the two evolutionary lineages have the same rates; 

• DNA sequences are at the compositional equilibrium when they start to diverge (nucleotide frequencies are 
constant). 

We will see that even with relaxing the last two hypotheses some calculations can be performed, but it is worth noting 
that compositional equilibrium, if the last assumption is verified, is maintained during the course of evolution. 
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Denoting with fi the compositional equilibrium frequency of the nucleotide i with i G {A,T, G,C} and with = 
Ti^j the substitution rate from nucleotide j to i in the unit time. The distance between two sequences, can now be 
defined as 



d = 2<^/,M. = 2^^/, ^ r,, . (1) 

Since 1969, when Jukes and Cantor proposed their first one-parameter model for nucleotide subsitution in DNA, 
many different models of increasing complexity have been published. The general 4-state Markov model has 12 
independent parameters, G12 in figUl (for a review see Zharkikh j^). This number, and consequently the model 
complexity, can be decreased by further conditions on the parameters, leading to a plethora of different models. A 
possible choice is to take into account the property of no strand-bias, explained in figQ] It was introduced by Sueoka 
in 1995 3] and we generally refer to it as type 1 parity rule or PRl. This rule is easily understood thinking that, 
scoring the substitution on one strand, the same substitution can be obtained in two ways: A ^ C is observed also if 
on the opposite strand T — > G. 

This means that we cannot discriminate substitutions between two bases from those between their complementary 
bases. In symbols: 

'''ij = Tljj (2) 

where the bar means complementary nucleotide: A = T and viceversa. And C = G similarly. 

The number of independent parameters is then halved, so that the following substitution rates can be introduced: 



a = 


rAT 


= rjA 


b = 


rAG 


= rjc 


c = 


rcT 


= rcA 


d = 


TAC 


= rjG 


e = 


rcA 


= rcT 


f ^ 


rcG 


= rcc 



The notation introduced here is consistent with the one previously used by Sueoka and Lobry '3| 
Equilibrium frequencies for such a model are easily derived from the master equations: 
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where qt denotes in general the probability of state i. 
These frequencies are given by: 



f _ „oc _ „oo _ 1 b + d 
Jl = - 9t - o 



h = Iq - Qc 



2b + c + d- 

1 c + e 
2b+c+d 



(3) 



The intrinsic symmetry of the model is evident. In this framework, in other words, there is only one independent 
frequency, the other being deduced by the normalization condition 2/i + 2/2 = 1. We now stress the fact that this 
is valid in a single strand {type 2 parity rule or PR2). If PRl is satisfied, then as a consequence the frequency of a 
nucleotide in a strand must be equal to that of its complement in the same strand. 

In the following we will resume some general results regarding PRl algebra showing that, in many cases, it is not 
possible to reconstruct the supposed underlying mutation pattern because the independent parameters outnumber 
the possible independent observations. 



II. MATERIALS AND METHODS 



In this section we will give some results regarding the model introduced above, focusing on the number of actual 
independent possible observations. 
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A. General model 

Given the substitution matrix R[4,4], whose entries are the mutation rates per nucleotide per unit of time, one can 
deduce the evolutionary matrix P[4_4](t), whose entries Pij{t) represent the probabiHty of finding at a certain site the 
base i at time t, given the base j at t = 0. Yet the divergence matrix X[4 4](t) can be deduced, whose entries Xij{t) 
are the mutual probability of finding at time t the base j in a sequence, given the base i at the same site of the other 
sequence. Obviously, if the substitution pattern is the same for both sequences, it results in Xij(t) — Xji{t). 

It is worth noting that the divergence matrix at initial time is nothing but the diagonal matrix with nucleotide 
frequencies on the diagonal. 

The result of an evolutive process can be synthetically represented as an initial diagonal divergence matrix, multi- 
plied on the left and on the right by a certain number of substitution matrices (corresponding to the generation steps 
in the two evolution lineages), producing a final matrix 

x(i) = R'^ . . . r^r; x(o) r1r'2---r; 

X{t) = P' X(0) P' (4) 

4 

X^ji*) = ^P\k{t)fkP]k{t) 
k=l 

where the substitution matrices can, in principle, all be different. 

The entries of the divergence matrix are the experimentally observable quantities. 
In our case the substitution matrix is R[4,4]: 
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1 — a — c — e 
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c 
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1 — a — e — c 
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d 


1-b-d-f 
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C 


d 


b 


f 


1-d-b-f 



obtained under the hypotheses of no-strand-bias, I.E. PRl. 

B. Non identifiability of some models 

In the following we show that the mathematical properties of the PRl algebra are such that, dealing with the 
general model, the parameters to estimate outnumber the possible independent observations, so that the model is 
untractable. As seen in eq.Q 

X{t) = P' X(0) P*. 

Now, several cases are possible, depending on whether P' P or not. In the following, we will assume that X(0) is 
already at compositional equilibrium, I.E. 

9a = 'Zt = /l = XAA{t = 0) - XTT{t = 0) 

(5) 

9c = 9g = h = Xccit = 0) = XGG{t = 0) 

1. p' = p 

As P' = P it is clear that X(i) is symmetric (X(t) — X'(t)). We have to estimate 6 parameters (6 mutation rates) 
and we have only 5 independent observations. This happens because of the symmetry Xij = Xji, the normalization 
conditions and because x^, = x^^. In more detail: 



XAG = XqA = XtC = XCT 
XAC = XCA = XtG = XgT 
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xat = xta 

XCG = XgC 

xaa = xtt 
xcc = xgg 



Where xaa = xtt and xcc = xcc can be deduced by the other four using the normahzation xij = fi) and the 
equihbrium frequencies. We find that xag, xac, xat, xcg and one equihbrium frequency are the only independent 
observable quantities. 



2. P' / P 



In this case mutation rates double becoming 12; so we have 12 parameters to calculate. Independent observations, 
on the other hand, increase up to 7, because of the lack of the symmetry Xij = Xji. Still the model is intractable. 



C. Reversible PRl model 



In this section we will deal with one of the previous models, the simplest one where P' = P. In this case simple 
calculations lead to an analytical expression for the divergence matrix, but the model remains intractable. Yet we will 
see that by the imposition of a certain property the model becomes tractable, and a way to estimate the parameters 
for a real data set will be proposed. 

In the following we will assume again that the initial divergence matrix is already at compositional equilibrium. 
Further, we will treat the evolutionary process as a continuous time process, being the time since the divergence very 
long. This allows us to write the following equations to solve the problem. The expression for the evolutionary matrix 
is 

P(t) = exp{Rt}; (6) 
as it is the solution of the differential equations (see Rodriguez et al. @) 



d9{t) 



dt 

^ = Ep..W^... (8) 

While the divergence matrix is given by 



P{t)R (7) 



k=l 



X{t) = P'{t)X{t^O)P^{t) (9) 

4 
k=l 

It is easily verified that, if P' = P, then Xij{t) = Xji{t). 

Now, the expressions for Xij{t) (the observables) can be inverted to obtain the rates and then the distance. 
The strategy could be: 

• solve the model, that is find the {t) as a function of rates; 

• invert the above equations to get an expression for the rates; 

• substitute the observed quantities Xij in order to have a numerical estimation of the rates; 

• use these estimates to obtain the distance. 

The expressions for can be deduced in a manner analogous to that proposed by Takahata & Kimura in 1981 
who deal with a slightly less general model than this (model TK5 in fig[21). In this way we get an expression for every 
entry of the divergence matrix, but with five independent expressions, as stated above. We repeat here the reasons: 



5 



• the symmetry of the matrix 

• the intrinsic symmetry of the model Xij = x^f, 

• the normaUzation conditions Xij — fi. 

Thus, we can write down the entire divergence matrix by means of the foUowing quantities: 

P ^ xag = xga = xtc = xcT 
R = xac = XCA = xtg = xgt 

Ql = XAT = Xta 
Q2 = XCG = XGC, 

51 = XAA = XtT 

52 = xcc = xgg 

Where, as stated above. Si and S2 can be deduced by the other four using the normahzation and the equihbrium 
frequencies. We find that P, i?, Qi, Q2 and one equihbrium frequency are the only independent observable quantities. 



D. Solution of the model 



Deriving an analytical expression for the divergence matrix is quite an easy task following g. Let's consider for 
example the element Xf\c] its derivative will be 

dx^c d{qAqc) . , . 

^-^^=9C9A + 9A-?C. (11) 

It is worth giving a brief explication for this. We said that we are considering the two lineages at compositional 
equilibrium at the initial time, so one would naturally say that qi = 0, and so the above equation. Stating that we 
are at compositional equilibrium means that sampling the whole considered sequence nucleotide frequencies fi 
don't change (apart from finite-size fiuctuations) . It does not mean that there is no mutation at all on each site; had 
this been the case, there would be no evolution to study. The probability for each nucleotide to mutate into another 
is given by the master equation, and this is why we can write Xij as qi times qj, take the derivative, and reexpress in 
terms of other qiqj products, I.E. other X entries. 
An example of derivative would be, for example, 

q/\ = [dqc + bqc + aqj) - (a + c + e)qA. (12) 

Substituting this and the analogue for qc in ea. Hll() and doing the same for all X entries we obtain a set of linear 
coupled first order differential equations which can be diagonalized and solved. 
More detail on the derivation is reported in the appendix El 



E. Reversibility 



Until now we have stated that it is possible to write the divergence matrix for this model, but it would be of no 
use because we could never invert five expressions and obtain six independent rates as functions of the matrix entries. 
What can be done is to reduce the number of independent parameters by adding a relation between them. Many 
choices are possible. One could be, following 0], a = /. Another possible choice is to make the model time reversible. 
We remember that time reversibility is satisfied when 

Ptjfj=Pjift (13) 

where pij are the entries of the evolutionary matrix and fi the equilibrium frequencies. It is possible to demonstrate 
that this property is equivalent to the detailed balance (see appendix^ which reads 



^ij^j ^jifi ^^iJ- 



(14) 
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In our model detailed balance holds if and only if 

be = cd. (15) 

This can be deduced by inspection of equilibrium frequencies expressions, or by a simpler rule 0, reported here in 
appendix El A general version of reversible model has been studied by Yang 3, who pointed out its ability of fitting 
the data better than other models. Gu and Li 9] have shown its robustness against violation of time reversibility. 

III. RESULTS AND DISCUSSION 

A. Estimation of the substitution rates 

Due to the complexity of the expressions coming from this model, it is hard to think that one can find an analytic 
way to invert them and express the rates as a function of the observables. Therefore we chose a statistic way to 
perform this inversion, based on the test. We write the ^ 

■ ■ •^^,7 ■ ■ ^2,7 

It is easily seen that this quantity is always non-negative, being zero when Xij = cc^j , I.E. when the model perfectly 
fits the observations. Clearly, by performing a minimization on it we look at the same time for the best parameters. 
In this contest trying to minimize the ^.s a function of six parameters would outcome in a complete failure, the 
algorithm would wander among the infinite number of equivalent solutions. Enforcing the reversibility makes the 
estimation possible, as it will be shown below. 

B. A realistic example 

As an application, we started from the multiple alignment of rRNA sequences used in [Tol |. The observed divergence 
matrix (unnormalized) between Xenopus and Homo is reported here below. 



Xenopus 
A T G CT 
A 647 1 17 2 
Homo T 3 523 11 18 
G 17 9 903 28 
C 8 21 25 691 



By changing parameter values over 6 magnitude orders we found that the criterion was well shaped with only 
one global minimum (fig. 13)1 . A systematic exploration of all possible pairs of parameters showed that there were 
no strong structural correlations between parameters, except between b and c (fig. As a consequence, parameter 
values are easily estimated using standard non linear minimizing tools (note that it is advisable to enforce parameter 
positivity during optimisation) . This example showed that parameter can be estimated in practice from a realistically 
sized dataset. 



C. Discussion 

The most general model of evolution at the DNA level has 12 parameters and this is too much for practical purposes. 
If we try to simplify it by enforcing some parameter to be equal, then the number of possible sub-models rapidly 
increases because many ways of doing it are possible. At the opposite side we find the only model which requires all 
the parameters to be equal (JC). 

It is clear that the number of published models in the literature doesn't cover all possible ones, and only those 
coming from some biological or mathematical justifications have been explored. 
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Under PRl hypothesis, we are dealing with no strand-bias models whose most general form has 6 parameters. We 
do not claim that models of this class are the best in any way, but that they are an interesting starting point. An 
important property of these models is their convergence towards PRE state even if substitution rates are modified 
during the course of evolution 111. PR2 state is a strong assumption and strand asymmetry has been observed in 
many cases. But, as PR2 is usually observed at a genome scale level 41, the hope is that, on average, with local 
deviations from PRl hypothesis canceling out, this class of model is not too bad an approximation. The biological 
motivation leading to the no strand-bias models has an important mathematical consequence, so, if it is biologically 
reasonable to study these models, one must be aware of the fact that the symmetry involved inexorably reduces the 
number of independent observations, making the model mathematically intractable. 

D. Conclusion 

As we have shown in section FlI B II comparing the number of unknowns to possible independent observations there 
is definitively no hope to estimate the 6 parameters of the general form of the no strand-bias model from pairwise 
DNA sequence comparisons. There is no unique solution to a system of M equations in A'^ > Af unknowns, in our 
case there is an infinite number of way to choose the six rates a, b, c, d, e, / in order to satisfy the five independent 
equations defining the matrix X. This result is extremely unpleasant because it corresponds to the most common 
situation with experimental data from present day DNA: fossil DNA data are scarce and from a relatively recent past. 
We clearly need further simplifications. 

We have exhibited here an example of a model, noted RNSB in figure 2, that combines the properties of reversible 
models and no strand-bias models. It is important to note that this model has still 5 parameters free because if the 
intersection between the reversible model class and the no strand-bias class were only -say- 3 parameter free models, 
there would not have been much flexibility left for further research. We do not claim that this new RNSB model is 
the best intersection between the two classes. We just claim that the RNSB model proves that it's possible to do so 
with 5 free parameters, so that there is no bottleneck here for further theoretical work on the parametric forms for 
this class of DNA substitution models. 
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APPENDIX A: DERIVATION OF THE DIVERGENCE MATRIX 

In order to obtain the expressions for the divergence matrix we define (following the notation introduced above) 

X± =25i±2Qi 

y± =252±2Q2 (Al) 
Z± = 4P ± 4R. 

These expressions reduce the problem to six first order ordinary coupled differential equations. This system is 
block-diagonal, can easily be inverted and its solution is: 

X+ = uj[uj + {1 - uj)e^"*] 

Y+ = (1 -tj)(l -w + we^«*) (A2) 

Z+ = 2uj{l - uj){l - e^°^) 

and 



+ [Cc^ + /32(l-c.)]e^=* + 
9 

+ [a2w + ?7(l-w)]e^=* + 
Z_ = ^{-2(^-7)[aw-/3(l-c^)]e^i* + 

r 

+ [a{5 - 7 + ffV - /3(<^ - 7 - 5)(1 - ^)]e^^* + 
+ [a(<5 - 7 - ffV - /3(<^ - 7 + 5)(1 - ^)]e^^*} 

where 



a = 


c — e 






P = 


b-d 






7 = 


2a + c + e 






5 = 


b + d + 2f 






UJ = 


2/1 = 2/^ = 


= 2fT 




Ao = 


-2(6 + c + 


d + e) 




Ai = 


-(2a + 6 + 


c + d + e ^ 


f 2/) 


A2 = 


Ai +g 






A3 = 


Ai - g 






a = 








C ^ 




- 7 + 5) 4 




■f] = 


lis-i)is 


-7-5)4 


- q;/3 



Combining all these, the entry for the divergence matrix are obtained. 

APPENDIX B: REVERSIBILITY AND DETAILED BALANCE 

We will show here the equivalence between time reversibility and detailed balance. 

1. DETAILED BALANCE => TIME REVERSIBILITY 

Let's just remind that 

P{t) = exp{Ri}, 

which can be developed as 

P{t) =I+Rt+^R2t2 + ..., 

or 

k 

Equation (|B2(I can be also written as: 

00 (n) 
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where 



felfc2---fcn-l 

for n > 2 



s-"^ = Tij, for n = 1. 



(B4) 



Now we will show that, if detailed balance is satisfied, then 

J J 



^^h^^h: V*,j,n. (B5) 



In fact, exploiting detailed balance, 

"^"■"/j = ri^kr---rk„_^,jfj (B6) 



ki---kn-i 



becomes 



and finally 



Sfel---fc„_l ^i.fcl ' ' ' '''j;k„-lfkn-l — 



»-fei,irfe,,fei •■•rj-fc„_J,. (B7) 

ki ■ ■■kn — i 



Reordering all the factors 



ki---kn-i 

ki •••kji — i 

As the sum is performed on indices ki ■ ■ ■ fc„_i the expression in IjBSp is equal to s"-^'^ fi for all n > 2. So we have ljB5|) 
for n > 1, and it is evident for n = \. Further, as 5ijfj — Sjifi, we obtain Pijfj ~ Pjifi Q- E. D. 

2. DETAILED BALANCE <= TIME REVERSIBILITY 

Let's rewrite the formula 

k 

Let's compute the time derivative oi pijfj] if time reversibility holds it will be equal to the time derivative of Pjifi- 
From the formula (|B9() . as equilibrium frequencies don't depend on time 

^(K. W/.) = = YP^k{t)r,,f,. (BIO) 



But 

dpij{t) 



dt 

k 



as P and R commute (evident from the solution). The second expression in IBlOp can be written as 
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Because of the time reversibility the last expression in (|B11|) becomes 



(B12) 



k k 



Finally 



4 fei = fi 



dt 



^P]k{t)rkih- 



(B13) 



k 



Subtracting the ljB13p from the ljB12p . which are equal, and keeping in evidence Pjk{t) we finally obtain 



^Pjk{t){rikfk - Tkifi) = 0, 



(B14) 



k 



and the detailed balance is satisfied Q. E. D. 



APPENDIX C: DETAILED BALANCE: SIMPLE CHECK 



A nice property of detailed balance is that there exists a very easy way to state if it holds, even without calculating 
equilibrium frequencies. Until now we have seen that the detailed balance is fulfilled when the equilibrium frequencies 
and the mutation rates (from which the former depend) cancel every term in the master equations. 

Another way to check the detailed balance is to consider three states in the system and the rates connecting 
them. If the product of the three rates which takes from a state to itself "clockwise" is equal to that calculated 
"counter-clockwise", then the detailed balance holds. If we have three states i, j, k then the above property will read 



[1] Kimura M., 1968. Evolutionary rate at the molecular level. Nature 217, 624-626. 

[2] Zharkikh A., 1994. Estimation of evolutionary distances between nucleotide sequences. J. of Mol. Evol. 39, 315-329. 

[3] Sueoka N., 1995. Intrastrand parity rules of DNA base composition and usage biases of synonymous codons. J. of Mol. 

Evol. 40, 318-325. Errata 42, 323 
[4] Lobry J. R., 1995. Properties of a general model of DNA evolution under no-strand bias conditions. J. of Mol. Evol. 40, 

326-330. Errata 41, 680. 

[5] Rodriguez F., Oliver J. L., Man'n A., Medina J. R., 1990. The general stochastic model of nucleotide substitution. J. of 
Theor. Btol. 142, 485-501. 

[6] Takahata N., Kimura M., 1981. A model of evolutionary base substitution and its application with special reference to 

rapid changes of pseudogenes. Genetics 98, 641-657. 
[7] Peliti L.; Appunti di meccanica statistica, Bollati Boringhieri (Torino, Italy, 2003) 
[8] Yang Z., 1994. Estimating the pattern of nucleotide substitution. J. of Mol. Evol. 39, 105-111. 

[9] Gu X., Li W.-H., 1994. A general additive distance with time-reversibility and rate variation among nucleotide sites. Proc. 

Natl. Acad. Sci. USA 93, 4671-4676. 
[10] Gouy M., Li W.-H., 1989. Phylogenetic analysis based on rRNA sequences supports the archaebacterial tree rather than 

the eocyte tree. Nature 339, 145-147. 
[11] Lobry J. R., Lobry C, 1999. Evolution of DNA base composition under no-strand-bias conditions when the substitution 

rates are not constant. Mol. Biol. Evol. 16, 719-723. 



^ik^kj^ji ^ij^jkTki- 



11 



c 

I 

A 



DNA 



I 

G 



G 

I 

T 



A 



no strand-bias condition means 



FIG. 1: Explication of the no strand-bias condition. If the rates for a certain substitution are the same on both strands of DNA, one can 
deduce the equivalence of this rate to the one between the complementary bases. 




FIG. 2: Hierarchy of DNA substitution models. Simplifications leading from a model to a simpler one are indicated by arrows. Only 
those directly referring to our discussion are drawn. This figure has been adapted from Robert Schmidt's work. 



Individual Influence of parameters near optimum 
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FIG. 3: shaped as a minimum over 6 orders of magnitude. 



Pairwise influence of parameters 
b and c near optimum 




0.009 0.010 0.011 0.012 0.013 0.014 



FIG. 4: Near the optimal values for the parameters, only 6 and c show a structural correlation. 



